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ABSTRACT 


The  linear  axisymmetric  triangular  element  is  used  to  study  the 
acoustics  of  axisymmetric  fluid  regions  subject  to  mixed  boundary  condtions. 
Asymmetric  acoustic  fields  are  represented  by  means  of  Fourier  Series 
expansions  in  the  circumferential  coordinate. 
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INTRODUCTION 


The  Finite  Element  Method  (Ref  l)  is  used  to  analyse  the  acoustic  field 
in  an  axisymmetric  fluid  region  subject  to  mixed  boundary  conditions.  The 
simple  ’linear'  axisymmetric  triangular  element  is  used  Li  the  analysis,  and 
asymmetric  acoustic  fields  are  represented  by  means  of  Fourier  Series 
expansions  in  the  circumferential  coordinate.  Although  this  note  is 
essential1/  written  from  an  acoustic  point  of  view,  the  necessary  ingredients 
are  present  to  facilitate  numerical  solution  of  a  variety  of  hydrodynamic 
potential  flow  problems,  eg,  non-axisymmetric  flow  over  bodies  of  revolution. 
Section  2  discusses  the  formulation  of  the  acoustic  problem  in  terms  of 
finite  elements  and  Sections  3  and  4  discuss  'free'  vibration  and  radiation 
problems.  Fortran  computer  programs,  described  elsewhere,  (Ref  2),  were 
written  to  solve  a  variety  of  problems;  some  examples  are  given  in  Section  5« 

2.  PROBLEM  FORMULATION 

Consider  the  acoustic  wave  equation  in  the  volume,  V,  enclosed  by  the 
surface,  Si 


div(grrl  <f>)  +^Q  -  \  ^  -  o  (l) 

where  <f>  is  the  acoustic  velocity  potential,  related  to  the  excess  pressure 
by  the  formula, 

p  =  p 

dt 

and  to  the  fluid  particle  velocity  by  the  formula, 

v  =  -  grad  4> 


C  is  the  sound  speed  in  the  fluid  and  p  is  the  density.  The  term,  Q, 
represents  sound  sources  which  may  be  active  in  the  fluid,  og,  a  point 
source  of  strength,  S,  say,  at  the  position  X  =  X^,  has  the  representation, 

Q  =  s  6  (X  -  xj 


The  boundary  conditions  are  assumed  to  be  of  nixed  type, 


dn 


+  <1  + 


(2) 


where  d  denotes  differentiation  along  the  outward  normal  to  the  fluid 

.  toB 

region. 

In  +he  usual  Finite  Element  sense,  'stiffness'  and  'mass'  matricc3  ere 


Preceding  page  blank 
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obtained  by  minimization  cf  the  following  functional  (Ref.  l)» 


•  ■if 


grad  grad  i  dv  -  /  (  Q  ”  ^  #Lv 

v  \  o  at8  / 


/  (  +  |  Ma  j 


Consider,  now,  the  axisymmetric  triangular  element  of  fluid}  a  cross 
section  is  shown  in  Fig.  1A,  referred  to  cylindrical  coordinates  (r,  0,  z) 
Let  the  variation  in  ^  over  this  element  be  given  by, 

<4  (r,  6,  z)  -  0n(r,  z)  cos  (n0)  (4) 
and  assume  that  (r,  z)  has  the  simple  linear  expansion 


4>n  (r,  z)  =  a1  +  a2  r  +  a^z 


It  is  not  difficult  to  show,  in  matrix  notation,  that 


+  (r,  0,  z)  =  [l  r  z]  [A-1]  ^  cos(n0) 


s™4  *  •  (to-  to  •  r  3e) 


cos(n0)  0  0  0  1  0  ^ 

»  0  oos(n0)  0  0  0  1  [A-1] 

0  0  sin(nQ)  -n  -n  -nz  k* 

r  r 

^ni,  are  the  values  of  (r,z)  at  the  nodal  points  (i,j,k) 

of  the  triangular  cross-section,  and 


[A-1]  =  {  F 


x  rA  'Vj 


zi  ‘  *k 


r*  “rJ 


Vi  “  ri*k 


*k  “  ai 


ri  “  rk 


X  =  +  r^z^-z^  +  r^z^Zj) 


Vj  “  rj*i 


2  a  Z 

i  i 


r 3  '  ri 


Substitution  of  (6)  and  (7)  into  (3)  and  then  minimizing  +  with 


4 


respect  to  ^  and  gives  the  equations, 


rdltA^J 


Q  2  Qn(r,s)oosn6 
q  =  qn(r,*)cosn0 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 


and  0  is  independent  of  6, 


r 


The  6-ir legration  (involving  sin*n0  or  cos*n8)  has  been  carried  cut 
in  (9)  through  (l3)j  A  denotes  the  cross-section  in  the  (r,z)  plane  and 
1  indicates  distance  alone  the  boundary  line  of  the  cross-section, 

is  equal  to  1  for  n  *  0,  and  en  “  2  when  n  -  0. 

The  assembly  of  individual  elements  to  form  a  region  (Pig  IB)  reflects 
the  fact  that  the  region  functional  is  the  sum  of  each  element  functional. 
The  assembly  is  essentially  an  organizational  problem  and  is  best  left  to  a 
computer  program.  A  matrix  equation  of  the  type. 


[M]  «L  +  [S]  +  [C] 


at 


I  j*ni  J  -  pniqj  -  j^niqj 


(15) 


results  where  is  a  column  vector  of  the  values  of  (r, z)  at  all 

structure  nodes  and  the  matrices  [C]  and  J  arise  entirely  from 

contributions  from  the  exterior  boundary  of  the  region. 

A  general  form  of  the  solution  of  (l)  subject  to  boundary  condition  (2) 
may  be  stated  in  the  following  manner.  Let  the  value  of  0  at  a  structure 
node  (i)  be  denoted  by  then  we  have  the  following  expansions  : 


<t>,  =  6  .  +  1.  .  cos(n8)  +  Z.  6  .  sin  (n0)  (l6a) 

ri  roi  n=1  rm  v  '  n=1  v  x  ' 

00  00  — 

Q(r,8,z)  »  Qo(r,z)  +  ^  (^(r.zJcosCne)  +  ^  Q^r^sinUe)  (l6b) 

00  00 

q(r,8,z)  =  qQ(r,z)  +  qn(r,z)cos(ne)  +  ng1  qn(r,z)sin(ne)  (l6c) 

P(r,0,z)  =  3(r,z)  independent  of  0.  (l6d) 


The  orthogonality  property  of  the  harmonic  functions  allows  reduction  from 
a  threi-dimensional  problem  (r,0,z)  to  a  series  of  two-dimensional 
problems  (r,z).  Assuming  that  the  Fourier  Series  (l6b)  and  (l6c)  may  be 
truncated  at  n  =  m,  then  2m  +  1  two-dimensional  problems  need  to  be 
solved  separately, 


1  involving  4>  . , 
01 

m  involving 
m  involving  <j>  . , 


where  i  runs  from  1  to  the 
number  of  nodes  in  the  region 
cross-section. 


The  total  potential  may  be  found  by  summation  of  the  series  (16a). 

It  should  be  noted  that  the  [ll]  and  [C]  matrices  are  independent  of  harmonic 
number,  and  the  [S]  matrix  for  the  is  identical  to  that  of 
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Part*  of  the  boundary  region  where,  q 
matically  set  to  zero  velocity,  [  d6 

\dn 

pressure  (0  -  0)  may  be  constrained  by  modification  of  the  ’assembled'  matrix 
equation.  The  integrals  involved  in  (9)  through  ( 1 3)  may  be  computed  from 
analytical  expressions  or  numerical  approximations;  the  problem  of  the 
singularity  at  r  ■  0  may  be  overcome  by  giving  the  relevant  nodal  r-coordinate 
a  small  positive  value  and  forcing  the  condition,  if  required,  of  0  -  0  at 
this  node. 

3.  PROBLEM  OF  FREE  YIBRATIOPS 


0  are  not  specified  axe  auto- 
and  parts  required  ;o  have  zero 


The  problem  of  free  vibrations  is  solved  by  setting  the  RHS  of  ( 1 5)  to 
zero  to  give  the  matrix  equation,  for  harmonic  number  n, 

£  -w8  [M]  +  [Si  +  =  0  (17) 

where  tine  variation  exp(-ikt)  has  been  assumed. 

The  relevant  boundary  condition  is,  from  (2), 

+  00  =  0  (IB) 

If  the  boundary  is  assumed  rigid  (or  pressure  release),  then  0,  and 
hence  [C],  vanishes  to  leave  (17)  as  a  standard  eigenvalue  problem. 

If  the  boundary  is  locally  reacting  with  acoustic  impedance, 

p/Vn  =  pc[  R(u)  -  iX(w)]  (l9) 

where  R  is  the  resistance  function  and  X  the  reactance  function,  then, 

0  =  K(X  -  iR)  /  (Xs  +  if),  K  =  tt/c  (20) 

and  the  eigenvalue  problem  is  non-standard.  However,  at  'low'  frequencies 
it  is  usually  possible  to  neglect  R,  and  the  ratio,  K/X,  is  independent 
of  frequency;  hence,  we  have  a  standard  eigenvalue  problem  for  many  cases 
of  practical  interest,  eg,  Helmholtz  resonators. 

4.  PROBLEM  OF  ACOUSTIC  RADIATIOII 

Consider  an  axisymmetric  surface,  S,  vibrating  with  known  normal 
velocity  (into  S),  Fig.  2A, 

Vs  =  Vn  (r,z)  cos(n9)  =  ”  (2l) 

s 

Surround  this  surface  by  a  sphere  of  radius,  R,  Fig.  2A,  such  that  R  is 
large  enough  to  allow  a  radiation  boundary  condition  on  the  sphere  surface, 
viz, 

^  ~  exp(iKR)  cos(n0)/R 

dd>  .  d±  -  (iX  -  1  /h)  4> 

dn  dR 
8 


on 


7 


(22) 

(23) 


“1 


Comparison  of  the  boundary  condition  (2)  with  (21 )  and  (23)  givest 


q  -  V^-,.;)  cos(n8)  on  the  vibrating  surface 
3  -  (t/R  -  iK)  on  the  large  sphere 


(24) 

(25) 


The  region  between  the  surface,  S,  and  the  sphere,  R,  is  divided 
into  triangular  elements,  Fig.  2B;  (24)  and  (25)  substituted  into  (9) 

through  (14)  then  give  the  eleme  nt  matrices  for  assembling  into  the  region 
matrix  ( 1 5 ) »  which  nay  be  solved  by  standard  techniques  to  yield  the 
potential  amplitude,  0^,  at  each  node.  The  above  analysis  is  repeated 

for  each  harmonic,  n,  to  give  the  general  potential  field, 


m  m  _ 

^  .  +  E  6  .  cos(n0)  +  Z  6  .  sin(ne) 
n=1  n=1 


(26) 


caused  by  boundary  vibration, 

m  m  _ 

Vc  =  V  (r,z)  +  E  V  (r,z)  cos(n0)  +  Z  V  (r,z)sin(n0)  (27) 
u  o  .  n  .  n 

n=1  n=1 

If  active  sound  sources,  Q,  are  present  in  the  fluid  the  procedure 
should  be  as  follows,  in  the  usual  acoustic  scattering  notation. 


Let  *  =  [^(r,z)  +  *s(r,z)]  cos  (n0)  (20) 

-.7here  ^i(r,z)cos(n6)  is  the  incident  potential  associated  with  the  active 

source  distribution,  Q,  and  ^g(r,z)cos(n9)  is  the  scattered  potential. 

Assuming  that  the  surface,  S,  is  passive  with  local  impedance  z(S),  we 
nave  the  follojing  boundary  condition  on  S, 


or. 


(4r  +  “A(s))  + 

k  +  iK  ^./z^A 

dn  \  dn  / 

e  'a  ' 


+  iK  ^g/z(s)  =  0 


(29a) 

(29b) 


s  -  s 

The  boundary  condition  to  be  applied  on  R  is  simply  (23), 

+  (l/R  -  iK)  =  0 
dn 


(30) 


The  problem  is  thus  formulated  in  terns  of  the  scattered  potential  <f>„, 
with  boundary  conditions, 


q  =  +  iK  ^/z(S) 

dn 

s 

3  =  iK/z(S) 


on  S 


(31) 
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1 


q  -  0 

0  -  (lA-  iK) 


on  R 


(32) 


The  source  tern,  Q,  is  not  needed  explicitly  in  the  element  properties  ( 1 2) , 
its  value  is  implied  in  the  boundary  condition  ( 31 )  • 


5.  SOME  NUMERICAL  EXAITLSS 

(a)  Consider  a  region  of  fluid  contained  in  the  space  between  two 
concentric  finite  cylinders.  A  cross-section  divided  into  triangular 
elements  is  shown  in  Fig.  3A;  there  are  00  elements,  the  bandwidth  of  the 
problem  is  13  and  the  total  number  of  degrees  of  freedom  (nodes)  is  55* 

The  boundaries  of  the  region  are  assumed  rigid.  The  natural  frequencies 
of  the  region  are  given  by  values  of  V7  where 

i  =  [A  Jn(Kr)  +  BYn(Kr)  ]  cos  ^  cos(n6)  (33) 

and  d£  is  to  vanish  at  r-a,  and  at  r=b. 
dr 

mV\ 

l2  ) 


n  identifies  the  circumferential  variation  of  $ 
m  identifies  the  axial  variation  of  0 
S  will  identify  the  radial  variation  of  0 


a 

=  8 

cms 

b  -  10  1-5 

c  -  150000 

n 

m 

S 

Finite  Element  Hz 

— 

Analytical  Hz 

* 

0 

1 

0 

15052 

15000 

0.3 

0 

2 

0 

30416 

30000 

1.4 

0 

0 

1 

3Q397 

37572 

2.2 

0 

1 

1 

41753 

40456 

3.2 

0 

3 

0 

46420 

45000 

3.2 

0 

2 

1 

50720 

48080 

5.5 

1 

0 

0 

2^59 

26>58 

0.0 

1 

1 

0 

15285 

15234 

0.3 

1 

2 

0 

30532 

301 18 

1.4 

1 

0 

1 

38481 

37667 

2.2 

1 

1 

1 

41843 

40544 

3.2 

1 

3 

0 

46495 

45078 

3.1 

1 

2 

1 

50791 

48154 

5.5 

2 

0 

0 

5315 

5315 

0.0 

2 

1 

0 

15963 

15914 

0.3 

2 

2 

0 

30877 

30467 

1.3 

2 

0 

1 

38761 

37964 

2.1 

2 

1 

1 

42099 

4O811 

3.2 

2 

3 

0 

46723 

45313 

3.1 

2 

2 

1 

51004 

_ mi2. _ 

,■■2-8 

The  above  Table  shows  good  agreement  between  the  Finite  Element  natural 
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frequencies  and  those  computed  from  the  a.  alytical  expression,  (33) I  good 
agreement  was  also  obtained  for  the  node  shapes.  The  fall  off  in  aocuracy 
as  m  and  S  increase  reflects  the  finite  division  of  the  fluid  region. 
Natural  frequencies  and  mode  shapes  were  also  obtained  for  the  above  geometry 
with  pressure  ~  release  boundary,  arid  for  rigid  and  pressure-releaso 
boundaries  of  a  finite  tube  of  fluid;  again  good  agreement  was  obtained 
between  Finite  Element  and  Analytical  results. 

(b)  Consider  the  Helmholtz  resonator  shown  in  Fig.  3B,  with  boundary 
assumed  rigid  except  where  the  opening  communicates  with  the  external  fluid. 
The  boundary  condition  at  the  opening  is  that  of  (20)  with  R<<X  and  set 
to  zero, 


3  -  K/X  (34) 

and  it  is  assumed  that  this  ratio  is  independent  of  frequency  in  the 
frequency  range  of  interest.  The  Helmholtz  resonance  corresponds  to  the 
first  natural  frequency  of  the  mode  n  *  0.  The  Finite  Element  natural 
frequency,  ascertained  in  cases  (i)  and  (ii)  below,  is  to  be  compared  with 
122  Hz  measured,  121  Hz  computed  from  an  ’improved  formula'  and  133  Hz  from 
the  classical  formula  which  is  seriously  in  error  for  this  shape  of 
resonator,  (Ref  3)  : 

(i)  If  we  choose  X  (  -  8bK/3w)  as  the  well  known  formula  for  the 

reactance  function  of  a  piston  vibrating  in  an  infinite  baffle  the 
Finite  Element  computation  gives  117  Hz. 

(ii)  If,  from  a  Finite  Element  Radiation  computer  program,  we  compute 
the  reactance  function,  X,  of  the  geometry  of  Fig.  3B  vibrating 
with  constant  velocity  at  the  opening  and  zero  velocity  elsewhere 

(X  -  1.34K),  we  obtain  a  Helmholtz  resonance  of  125  Hz. 

Bearing  in  mind  that  the  Finite  Element  method  gives  an  upper  bound  on  the 
exact  natural  frequency,  good  agreement  obtains  between  measured  and 
computed  frequencies. 

(c)  The  numerical  computation  of  acoustic  radiation  is  complicated 
by  the  necessity  of  choosing  the  radius,  R,  of  the  large  sphere.  If  the 
chosen  radius  is  too  small,  then  incorrect  results  will  be  obtained  due  to 
the  enforcement  of  a  radiation  boundary  condition  where  such  c  condition 
does  not  apply;  if  the  radius  is  too  large  then  the  computational  effort 
needed  for  solution  may  be  excessive.  In  addition,  detailed  attention  must 
be  given  to  the  results,  eg,  it  may  be  necessary  to  refine  the  triangular 
subdivision  where  rapid  changes  in  acoustic  pressure  occur,  and  in  many 
cases  it  is  advisable  to  plot  acoustic  pressures  at  centroidal  rather  than 
nodal  points. 

Consider  the  geometry  of  Fig.  2A,  a  cylinder  of  length  20  cm  and 
radius  5  cm  with  stationary  hemispherical  end-caps.  The  region  between  the 
surface  and  a  sphere  of  radius  3B  cm  was  divided  into  392  elements  involving 
225  nodal  points  and  a  bandwidth  of  33;  on  the  surface  15  nodal  points 
were  chosen.  Fig  4  shows  nodal  acoustic  pressures  obtained  when  the 
cylinder  surface  pulsates  (n-0)  at  10,000  Hz  with  velocity  amplitude 
0.001  cm/sec  into  a  fluid  of  density  1.0  and  sound  velocity  150,000  ca/seo. 
Fig.  5  shows  acoustic  pressures  at  0  =  0  when  the  cylinder  rigidly  oscillates 
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(n-l).  For  comparison,  results  obtained  by  means  of  fitting  spheroidal 
wavefunctions  (Ref  4)  are  given.  Important  acoustic  parameters  are  the 


radiation  resistance  function, 
as  follows, 

R, 

and 

the  reactance  function,  X, 

defined 

R*al  ^  j 

f  pV*ds  - 
S 

1 

2 

Rp  J 

f  W*ds 

S 

(35) 

Imag.  ^  j 

f  pV  ds  - 

1 

2 

XpC 

f  W*ds 

J 

S 

4 

where  *  denotes  complex  conjugate.  The  first  of  (35)  is  the  acoustic 
power  radiated  and  the  second  is  the  'wattless'  power. 


Harmonic 

Finite 

Element 

Spheroidal 

Wavefunctions 

Number 

R 

X 

R 

X 

0 

0.911 

0.284 

0.939 

0.274 

1 

0.951 

0.350 

0.987 

0.352 

Figs.  4  an.i  5  and  the  shove  table  show  reasonable  agreement  between  the 
Finite  Element  method  and  the  series  expansion  method  using  spheroidal 
wavefunctione. 

6.  SUMMARY 

Tha  Finite  Element  method  has  been  applied  to  the  acoustics  of 
axisymmetric  fluid  regions  subject  to  mixed  boundary  conditions; 
asymmetric  fields  being  represented  by  means  of  Fourier  Series  expansions. 
Numerical  results  suggest  that  good  accuracy  may  be  obtained  for  interior 
problems  with  little  judgement  from  the  user  of  the  computer  programs. 
Acoustic/potential  flow  problems  involving  circular  ducts  of  varying 
cross-section  and  longitudinal  impedance  can  easily  be  solved  by  the 
method  presented.  The  exterior  (radiation)  problem,  however,  requires 
detailed  study  in  order  to  obtain  correct  sound  levels,  especially  as 
frequency  increases.  It  may  be  possible  to  overcome  this  problem  by 
using  curved  elements  to  represent  boundaries  and  simple  wavefunctions 
to  describe  the  potential  field  of  an  element.  However,  radiation  problems 
intractable  by  conventional  standards  may  be  solved  with  good  accuracy  if 
close  attention  is  paid  to  the  triangularization  of  the  region  concerned, 
eg  radiation  from  the  opening  of  a  circular  duct  and  radiation  from  a 
finite  ring  of  arbitrary  cross-section. 
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FIG.  2A  CROSS  SECTION  OF  VIBRATING  SURFACE  SURROUNDED 
BY  SPHERE  OF  RADIUS  R 
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FIG.  2B  THE  REGION  OF  FIG.  2A  DIVIDED 
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•  FINITE  ELEMENT 
—  SPHEROIDAL  WAVE  FUNCTIONS 


FIG.  4  PRESSURE  AMPLITUDE  AT  10000H?.  N  =  0 
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SPHEROIOAL  WAVE  FUNCTIONS 


FIG. 5  PRESSURE  AMPLITUDE  AT  10000Hz,  N  =  1 
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